function conditionNumber = calculateConditionNumber( H,precision,ErrTol, MaxIts, verbose,k)
%CALCULATECONDITIONNUMBER Summary of this function goes here
%   Detailed explanation goes here
    n=(size(H.filter,1)+1)/2;
    dim=ndims(H.filter);

    if nargin < 6
		%use strict default parameters
		k=100; % ??? calculateConsecutiveMatrixApplication(size(H.filter,1));
		%fprintf('use strict default parameters for calculateConsecutiveMatrixApplication, k = %d\n',k);
    end

    L=3*n-2;
    
    if dim == 2
        Hhat=fftn(H.filter,[L L]);
        W=randn(n,n);% W is distributed N(0,1) iid
        %W=1./(2*(rand(n,n)-0.5)); W=W-sign(W); %U[0,1] => N(0,1)
    elseif dim == 3
        Hhat=fftn(H.filter,[L L L]);
        W=randn(n,n,n);% W is distributed N(0,1) iid
        %W=1./(2*(rand(n,n,n)-0.5)); W=W-sign(W); %U[0,1] => N(0,1)
    else
        error('dim must be 2 or 3');
    end

    W=W/norm(W(:));% and normalized
    
    v=W;
    rho=1;
    for j=1:k
        if dim == 2
            v=feval(@applyHFilter2D,v,Hhat);
        elseif dim == 3
            v=feval(@applyHFilter3D,v,Hhat);
        else
            error('dim must be 2 or 3');
        end
        d=norm(v(:));
        rho=sqrt(d);
        v=v/d;
    end
    
    normH=rho;
    
    if nargin<5
       verbose = 0;
    end

    if nargin<4
       MaxIts = 20;
    end

    if nargin<3
       ErrTol = 1e-12;
    end
    
    v=W;
    rho=1;
    for j=1:k
        if dim == 2
            [v,flag,residual,iter] = ConjugateGradient('applyHFilter2D',v,{Hhat},ErrTol,MaxIts,zeros(size(v),precision),verbose);
        elseif dim == 3
            [v,flag,residual,iter] = ConjugateGradient('applyHFilter3D',v,{Hhat},ErrTol,MaxIts,zeros(size(v),precision),verbose);
        else
            error('dim must be 2 or 3');
        end
        d=norm(v(:));
        rho=sqrt(d);
        v=v/d;
        if flag
           %fprintf('Inversion did not converge. Residual error %-2.5e, iter = %d\n',residual,iter);
           %fprintf('iter=%d\trho=%e\n',j,rho);
        end
    end
    normInverseH=rho;
    conditionNumber=normH*normInverseH;